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analyze various observables for stable hadrons and hadron resonances produced in ultra-relativistic 

\Q '. 

. heavy ion collisions. Particle multiplicities are determined based on the concept of chemical freeze- 

out. Particles can be generated on the chemical or thermal freeze-out hypersurface represented by a 
parameterization or a numerical solution of relativistic hydrodynamics with given initial conditions 
and equation of state. Besides standard space-like sectors associated with the volume decay, the 



We have developed a fast Monte Carlo procedure of hadron generation allowing one to study and 



hypersurface may also include non-space-like sectors related to the emission from the surface of 



expanding system. For comparison with other models and experimental data we demonstrate 
the results based on the standard parameterizations of the hadron freeze-out hypersurface and 
flow velocity profile under the assumption of a common chemical and thermal freeze-out. The 
C++ generator code is written under the ROOT framework and is available for public use at 



http:/ / uhkm.jinr.ru/ 
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I. INTRODUCTION 



Ongoing and planned experimental studies of relativistic heavy ion collisions in a wide 
range of beam energies require a development of new event generators and improvement of 
the existing ones Especially for Large Hadron Collider (LHC) experiments, because of 
very high hadron multiplicities, one needs fairly fast Monte-Carlo (MC) generators for event 
simulation. 

A successful but oversimplified attempt of creating a fast hadron generator motivated by 
hydrodynamics was done in Ref. [2, 0, 0, 0]. The present work is an extension of this ap- 
proach. We formulate a fast MC procedure to generate hadron multiplicities, four-momenta 
and four-coordinates for any kind of freeze-out hypersurface. Decays of hadronic resonances 
are taken into account. We consider hadrons consisting of light u,d and s quarks only, but the 
extension to heavier quarks is possible. The generator code is written in the object-oriented 
C++ language under the ROOT framework p. 

In this article we discuss only central collisions of nuclei using the Bjorken-like and 
Hubble-like freeze-out parameterizations used in so-called "blast wave" [tJ and "Cracow" 
models 8(, respectively. The same parametrizations have been used in the hadron generator 
referred as THERMINATOR P| that appears however less efficient than our generator (see 
sections ITTllVll). 

The paper is now organized as follows. Sections HTHVl are devoted to the description of 
the physical framework of the model. In section IVI[ the Monte Carlo simulation procedure 
is formulated. The validation of this procedure is presented in section IV11I In section IV1111 
the example calculations are compared with the Relativistic Heavy Ion Collider (RHIC) 
experimental data. We summarize and conclude in section HXl 



II. HADRON MULTIPLICITIES 



We give here the basic formulae for the calculation of particle multiplicities. We consider 
the hadronic matter created in heavy-ion collisions as a hydrodynamically expanding fireball 
with the equation of state of an ideal hadron gas. 

The mean number Ni of particles species i crossing the space-like freeze-out hypersurface 
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a(x) in Minkowski space can be computed as 







Ni= / d 3 a,{x)j^{x). (1) 

Ja(x) 

Here the four-vector <i 3 <T M (x) = n fl (x)d 3 a(x) is the element of the freeze-out hypersurface 
directed along the hypersurface normal unit four- vector n^(x) with a positively defined 
zero component (n°(x) > 0) and d 3 a(x) = \d 3 a fM d 3 a fJj \ 1 ^ 2 is the invariant measure of this 
element. The normal to the space-like hypersurface is time-like, i.e. n^n^ = 1; generally, for 
hypersurfaces including non-space-like sectors, the normal can also be a space-like so then 
n^n^ = — 1. The four- vector jf{x) is the current of particle species i determined as: 

/ ^fi(x,p), (2) 

where fi(x,p) is the Lorentz invariant distribution function of particle freeze-out four- 
coordinate x = {x°,x} and four-momentum p = {p°,p}. In the case of local equilibrium 

f<{*,p) = fr( P -u(x);T(x)^(x)) = ^ e ^ {[v . u{x) J; i{x) y T{x))±V (3) 
where p ■ u = p^u^, gi = 2Jj + 1 is the spin degeneracy factor, T(x) and fii(x) are the 
local temperature and chemical potential respectively, u(x) = j{l,v} is the local collective 
four-velocity, = 1. The signs ± in the denominator account for the 

proper quantum statistics of a fermion or a boson, respectively. 
The Lorentz scalar local particle density is defined as: 

/d 3 p 
—p^(x)fi(x,p). (4) 

For a system in local thermal equilibrium, the particle density in the fluid element rest frame, 
where w* M = {1,0,0,0}, is solely determined by the local temperature T(x*) and chemical 
potential Hi{x*) for each particle species v. 

pf{T{x*\^)) = u^OO = / d 3 p*fr(f°;T(x*),p t (x*)); (5) 

the four-vectors in fluid element rest frames are denoted by star. 

In the case of local equilibrium, the particle current is proportional to the fluid element 
four-velocity: j^ix) = p° q (T(x), fii{x))u^ {x) . So the mean number of particles of species i 
is expressed directly through the equilibrated density: 

Ni= f d 3 a ll (x)u^(x)p^(T(x),p i (x)). (6) 

Ja{x) 



In the case of constant temperature and chemical potential, T{x) = T and fa(x) = fa, 
one has 

Ni = p?(T, fa) I <Pa,{x)u»{x) = p?(T, fa)V eS , (7) 

Ja(x) 

i.e. the total yield of particle species i is determined by the freeze-out temperature T, 
chemical potential fa and by the total co-moving volume V e s, so called effective volume of 
particle production which is a functional of the field of collective velocities u^ix) on the 
hypersurface a(x). The effective volume absorbs the collective velocity profile and the form 
of hypersurface and cancels out in all particle number ratios. Therefore, the particle number 
ratios do not depend on the freeze-out details as long as the local thermodynamic parameters 
are independent of x. The concept of the effective volume and factorization property similar 
to Eq. has been considered first in Ref. [ill, repeatedly used for the analysis of particle 
number ratios ( see, e.g., Ref. |l2j]) and recently generalized for a study of the averaged phase 



space densities 



window 
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13, 
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and entropy l^j]. One can apply this concept also in a limited rapidity 



The concept of the effective volume can be applied to calculate the hadronic composition 



at both chemical and thermal freeze-outs 



12[. At the former one, which happens soon after 



hadronization, the chemically equilibrated hadronic composition is assumed to be established 
and frozen in further evolution. The chemical potential fa for any particle species i at the 
chemical freeze-out is entirely determined by chemical potentials \x q per a unit charge, i.e. 
per unit baryon number B, strangeness S, electric charge Q, charm C, etc. It can be 
expressed scalar product: 

fa = Qifa (8) 

where qi = {B i: Si, Qi, Ci, ...} and Jl = {]2b, Ps, Pq, Pc, ••}• Assuming constant temperature 
and chemical potentials on the chemical freeze-out hypersurface, the total quantum numbers 
q = {B, S,Q,C, ...} of the selected thermal part of produced hadronic system (e.g., in a 
rapidity interval near y = 0) with corresponding V e R can be calculated as q = V e s Y^iPi^Qi- 
For example: 

n 

B = V cS J2pT(T,fa)B i , (9) 
i=i 



S = V cff PT{ T 1 Pi)Si 



(10) 



1=1 
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Q = V cS J2ti q (T, f i t )Q l . 



(11) 



The potentials Jl q are not independent. Thus, taking into account baryon, strangeness 
and electrical charges only and fixing the total strangeness S and the total electric charge 
Q, Jig and JIq can be expressed through baryonic potential Jib using Eqs. (p"0|) and fTTf . 
Therefore the mean numbers of each particle and resonance species at chemical freeze-out 
are determined solely by the temperature T and baryonic chemical potential JIb- 

In practical calculations, we use the phenomenological observation |l5j that particle yields 
in central Au+Au or Pb+Pb collisions in a wide center-of-mass energy range y/s NN = 
2.2 — 200 GeV can be described within the thermal statistical approach using the following 
parametrizations of the temperature and baryon chemical potential [poT ]: 



a = 0.166 ± 0.002 GeV, b = 0.139 ± 0.016 GeV" 1 , c = 0.053 ± 0.021 GeV~ 3 and d = 
1.308 ± 0.028 GeV, e = 0.273 ± 0.008 GeV" 1 . 



The particle densities at the chemical freeze-out stage are too high (see, e.g., |l2|) to 
consider particles as free streaming and associate this stage with the thermal freeze-out 
one. The mean particle numbers N} h at thermal freeze-out can be determined using the 

n 

following procedure j!2| . First, the temperature and chemical potentials at chemical freeze- 
out have to be fitted from the ratios of the numbers of (quasi) stable particles. The fitting 
procedure should account for the decays of all resonances as well as unstable particles in 
given experimental conditions (feed-down). The common factor, V^, and, thus, the absolute 
particle and resonance numbers can be fixed, e.g., from pion multiplicities. Within the 
concept of chemically frozen evolution these numbers are assumed to be conserved except 
for corrections due to decay of some part of short-lived resonances that can be estimated 
from the assumed chemical to thermal freeze-out evolution time. Then one can calculate the 
mean numbers of different particles and resonances reaching a (common) thermal freeze-out 
hypersurface. At a given thermal freeze-out temperature T t h these mean numbers can be 
expressed through the thermal effective volume V*g and the chemical potentials for each 



T((jl b ) =a- b[i 2 B - cf/ B , 



(12) 




(13) 
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particle species jjf 1 which can no more be expressed in the form /ij = q{jx valid only for 
chemically equilibrated systems. For a given parametrization of the thermal freeze-out 
hypersurface, the thermal effective volume (and thus all /4 h ) can be fixed with the help 
of pion interferometry data. 

In practical calculations we determine all macroscopic characteristics of a particle system 
with the temperature T and chemical potentials /i, via a set of equilibrium distribution 
functions in the fluid element rest frame: 

/fV°;7>;) = (27r)3exp( ^o_\ ]/r)±1 - ( 14 ) 

Eq. © for the particle number density then reduces to 

poo 

p?{T, in) = 4tt / dp*p* 2 fr(p*°; T, ta). (15) 
Jo 



Using the expansion 



/r(p* l T, ft) = ^ £(T) i+1 eMk^j^), (16) 



the density can be represented in a form of a fast converging series 



P?{T^i) = |^f;^exp(^ 2 (^), (17) 

71 k=l 

where Ki is the modified Bessel function of the second order. 

We assume that the calculated mean particle numbers iV, = pfKs correspond to a grand 
canonical ensemble. The probability that the ensemble consists of N{ particles is thus given 
by Poisson distribution: 

- (N-) Ni 

P(N t )=e W (-N t ) [ -^ r . (18) 

III. HADRON MOMENTUM DISTRIBUTIONS 

We suppose that a hydrodynamic expansion of the fireball ends by a sudden system 
breakup at given temperature and chemical potentials. In this case the momentum distribu- 
tion of the produced hadrons keeps the thermal character of the equilibrium distribution 



Similar to 
formula Il6l 



)qs. (JH), (j2J), this distribution is then calculated according to the Cooper- Frye 
P°^= / ^{xWfFiP'vWT,*). (19) 

d 6 P Ja(x) 



The integral in Eq. (JT§j) can be calculated with the help of the invariant weight 

d 6 N- 

W ^P)^P° 1 ^=n,(x)p^r(p-u(xy,T^ l ). (20) 
It is convenient to transform the four- vectors into the fluid element rest frame, e.g., 



(21) 



n* = n — 7(1+7) 1 (ri*° + n°)v 
and calculate the weight in this frame: 

W ff ,i(x,p) = W: 4 (x*,p*) = nl{x) V ^n\f°-T^i). (22) 

Particulary, in the case when the normal four- vector n^(x) coincides with the fluid element 
flow velocity u»(x), i.e. n*^ = u*^ = {1, 0, 0, 0}, the weight W^{x*,p*) = p*°f? q {p*°; T, /1) is 
independent of x and isotropic in the three-momentum p *. A simple and 100% efficient sim- 
ulation procedure can then be realized in this frame and the four-momenta of the generated 
particles transformed back to the fireball rest frame using the velocity field v(x): 

P° = l{p°* + vp*), ^ 
p = p* + 7(1 + 7)" 1 (p°* +p°)v. 

There are two well-known examples of the models giving n M (x) = u^(x): the Bjorken model 
with hypersurface tb = (t 2 — z 2 ) 1 ! 2 = const and absent transverse flow and the model with 
hypersurface th = (t 2 — x 2 — y 2 — z 2 ) 1 / 2 = const and spherically symmetric Hubble flow. In 
general case n^(x) may differ from u^{x) and one should account for the x — p correlation 
and the corresponding anisotropy due to the factor n M p M even in the fluid element rest frame 



IV. GENERALIZATION OF THE COOPER-FRYE PRESCRIPTION 

It is well known that the Cooper-Frye freeze-out prescription in Eq. ()19j) is not valid 
for the part of the freeze-out hypersurface characterized by a space-like normal four-vector 
n^. In this case \n°\ < \n\ and so p M n M < for some particle momenta thus leading to 
negative contributions to particle numbers. Usually, the negative contributions are simply 
rejected This procedure however violates the continuity condition of the flow piu^n^ 
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through the freeze-out hypersurface. Taking into account the continuity of the particle flow, 
the generalization of Eq. (fTUj) has the form ^| : 



p 



d 3 p 



a{x) 



d^(x)^(x,p)f^(T(x),tii(x)), 



(24) 



where 



(25) 



n^(x,p) = pi* 0(1 — \X(x,p)\) + u^(x) p ■ u(x) 0(\\(x,p)\ — I) 
X(x,p) = 1 — p ■ n(x) [p ■ u(x) n(x) ■ n(.x)] _1 , 
0(x) = 1 for x > 0, 0(x) = for x < 0. 

Passing to the fluid element rest frames at each point x and using Lorentz transformation 
properties of the quantities in Eq. (J24j), one arrives at the same form of the four- vector of 
particle flow as in the case of the freeze-out hypersurface with the time-like normal n M (x): 

d 3 p 



f(x) 



Po 



-^(x,p)fr(T(x),^(x)) = p?(T(x),H(x))vr(x) 



(26) 



Therefore the factorization of the freeze-out details in the effective volume in the case 
of constant temperature and chemical potentials, i.e. Eq. (JI|), is valid for any type of 
hypersurface |l3j . It follows from Eqs. ()24|) . (J2*5j) that the invariant weight in the fluid 
element rest frame has then the form: 



w;ax*, p * 



P*»< 1- 



p n 



+ p*°n*° 6 



p n 



- 1 



/fV°;T, Mi ). (27) 



For the time-like normal n M (x), Eq. ()27|) reduces to Eq. 

It is worth noting that though the bulk of particles is likely associated with the volume 
decay, the particle emission from the surface of expanding system, or formally, from a non- 
space-like part of the freeze-out hypersurface enclosed in Minkowski space, is essential for a 
description of hadronic spectra and like pion correlations at relatively large px j3- 



V. FREEZE-OUT SURFACE PARAMETERIZATIONS 

In principle, one can specify the fireball initial conditions (e.g., Landau- or Bjorken-like) 
and equation of state to follow the fireball dynamic evolution until the freeze-out stage with 
the help of relativistic hydrodynamics. The corresponding freeze-out four-coordinates x^, 
the hypersurface normal four- vectors n^(x) and the collective flow four- velocities u^(x) can 
then be used to calculate particle spectra according to generalized Cooper- Frye prescription. 



S 



This possibility is forseen as an option in our MC generator. In this paper, we however do 
not consider the fireball evolution, we demonstrate our fast MC procedure utilizing the 
simple and frequently used parametrizations of the freeze-out. 

At relativistic energies, due to dominant longitudinal motion, it is convenient to substitute 
the Cartesian coordinates t, z by the Bjorken ones 

r=(t 2 -z 2 ) 1 / 2 , rj = -ln t -±^ (28) 

and introduce the the radial vector r = {x,y} = {r cos0, r sin0}, i.e.: 

x^ = {t cosh??, r, t sinh rj} = {r cosh rj, r cos 0, r sin 0, rsinh^}. (29) 

Similarly, it is convenient to parameterize the fluid flow four-velocity u^(x) = 
7(x){l, v(x)} = 7(x){l, v r (x) , v z (x) } at a point x in terms of the longitudinal (z) and trans- 
verse (r) fluid flow rapidities 

1 l + v z (x) 1 1 + v r (x) coshr] u (x) 

Vu{x) = ~ In rT , pu{x) = - In — — — , 30 

2 1 — v z [x) 2 1 — v r (x) cosh rj u {x) 

where v r = \v r \ is the magnitude of the transverse component of the flow three- velocity 
v = {v r cos <f> u , v r sin <f) u , v z }, i.e. 

u^(x) = {cosh p u cosh r] u , sinh p u cos U , sinh p u sin U , cosh p u sinh 
= {(1 +u 2 r ) l/2 coshr] u ,u r , (1 + w 2 ) 1/2 sinh 7^}, 

u r = ■jVr = jrCOshrjuVr, 7 r = coshp^. For the considered central collisions of symmetric 
nuclei, U = 0. Representing the freeze-out hypersurface by the equation r = r(r), r, 0), the 
hypersurface element in terms of the coordinates rj, r, becomes 

,n dx a dx l3 dx' r , . 

d °» = 6 ^ drjdrdcf, dT]drd ^ (32) 

where e MQ!( g 7 is the completely antisymmetric Levy-Civita tensor in four dimensions with 
e 0123 = —eoi23 = 1- Particulary, for azimuthaly symmetric hypersurface r = t(t], r), Eq. 



yields 



121: 



d 3 a„ = r(r, ri)d 2 rd'n{ — d ^-- sinh r\ + cosh ri, — — cos 0, — — sin 0, ^— cosh 11 — sinh ri}. (33) 

t drj dr dr t dr] 

Generally, the freeze-out hypersurface is represented by a set of equations r = Tj(rj, r, 0) and 
Eq. (|32Jl should be substituted by the sum of the corresponding hypersurface elements. 



To simplify the situation, besides the azimuthal symmetry, we further assume the longitu- 



dinal boost invariance 2l( . The local quantities (such as particle density) are then functions 
of t and r only. The hypersurface then takes the form r = r(r), the flow rapidities r] u = n 
(i.e. v z = z/t), p u = p u (r) and Eq. yields 

cPcTfj, = t [r)d 2 rdn {cosh n, — ^ cos 0, — 4^ sin 0, — sinh?]}, 

d 3 (7 = |1 - (^i)2|i /2 r(r)rf 2 f^, (34) 
n^(x) = |1 — (^) 2 | _1 / 2 {cosh?7, ^ cos0, ^ sin 0, smli^}. 

Note that the normal four-vector n^ becomes space-like (n^n^ = —1) for \dr/dr\ > 1. 
For the simplest freeze-out hypersurface r = const one has 

d 3 a = Td 2 rdn, 

(35) 

n M (x) = {cosh r], 0, 0, sinhr?}. 

In this case the normal n^(x) is time-like (n^n^ = 1) but generally different from the flow 
four- velocity u^(x) except for the case of absent transverse flow (i.e. p u = 0). Assuming 
U = and the linear transverse flow rapidity profile (effectively taking into account a 
positive flow - radius correlation up to the radii close to the fireball boundary as indicated 
by numerical solutions of (3+l)-dimensional relativistic hydrodynamics, see, e.g., 2^|): 

Pu = (36) 

where R is the fireball transverse radius, the total effective volume for particle production 
at t = const is 

Kff = / d 3 a fl (x)u fl (x) = t j r rdr I d(f) drj = 

J <t(x) Jo Jo J »7 m i n 

= 2nrAn ( JL) (p^sinhp^ - coah/C* + 1), (37) 



1 values of the maximal transverse flow rapidity p™' 



where An = n max - n min . For sma 

Eq. reduces to V eS = nrR 2 An [12 1. 

We shall refer the above choice of the freeze-out hyper-surface and the flow four-velocity 

profile as the Bjorken-like parametrization or Bjorken model scenario for particle freeze-out 

with transverse flows 12111 . 

1 5 n 

We also consider so called Cracow model scenario |8j corresponding to the Hubble-like 
freeze-out hypersurface th = (t 2 — x 2 — y 2 — z 2 ) 1 ^ 2 = const and flow four- velocity 

u »(x) = x^/t h . (38) 
10 



Introducing the longitudinal space-time rapidity_n according to Eq. (|28|) and the transverse 
space-time rapidity p = sinh -1 (r/r#), one has |23j 

x M = r# {cosh rj cosh p, sinh p cos 0, sinh p sin 0, sinh 77 cosh p}, (39) 

t h = r B /coshp. Representing the freeze-out hypersurface by the equation t h = 
Th(t], p, 0) = const, one finds from Eq. 



(40) 



d 3 a = Tjj sinh p cosh pdndpdcf) = THdnd 2 f, 
rz M (x) = 

The effective volume corresponding to r = r H sinh p < R and 7? min < 77 < 77 max is 

r /" 2?r flmax 

Kff = / d 3 a fl (x)u^ l (x) —t h rdr / d0 / = irr H R 2 Ar]. (41) 

VI. HADRON GENERATION PROCEDURE 

Our MC procedure to generate the freeze-out hadron multiplicities, four-momenta and 
four-coordinates is the following: 

1. First, the parameters of the chosen freeze-out model are initialized. Particularly, 
for the models with constant freeze-out temperature T and chemical potentials pi, 
the phenomenological formulae ()12j) . ()13|) are implemented as an option allowing to 
calculate T and pi at the chemical freeze-out in central Au + Au or Pb + Pb collisions 
specifying only the center-of-mass energy y/s NN . In the scenario with the thermal 
freeze-out occurring at a temperature T th < T ch , the chemical potentials pf 1 are no 
more given by Eq. (jHJ). At given thermal freeze-out temperature T th and effective 
volume V^j , they are set according to the procedure described in section |TT| So far, 
only the stable particles and resonances consisting of u, d, s quarks are incorporated 
in the model. They are taken from the ROOT particle data table jfJS]. 

2. Next, the effective volume corresponding to a given freeze-out model is determined, 
e.g., according to Eq. (pT7j) or (jHJ) and particle number densities are calculated with 
the help of Eq. (|17jl. The mean multiplicity of each particle species is then calculated 
according to Eq. Q. A more general option to calculate the mean multiplicities, e.g., 
in the case of the freeze-out hypersurface obtained from relativistic hydrodynamics, is 

11 



the direct integration of Eq. (|24|). The multiplicity corresponding to the mean one is 
simulated according to Poisson distribution in Eq. (JTSJ). 

3. The particle freeze-out four-coordinates = {r cosh 77, r cos0, r sin0, r sinh?]} in the 
fireball rest frame are then simulated on each hypersurface segment r = Tj(r) according 
to the element d 3 a fl u fl = cPa^ = n* Q {r) |1 — (dr / 1 dr) 2 \ l l 2 T{r)d 2 rdn, assuming n* Q and r 
functions of r (i.e. independent of n,<p), by sampling uniformly distributed n in the 
interval [^ m in, ^max], m the interval [0, 2tt] and generating r in the interval [0,i?]) 
using a 100% efficient procedure similar to ROOT routine GetRandomQ. In the 
Bjorken- and Hubble-like models: r(r) = tb = const, Uq = coshp u = 7 r and |1 — 
(dr 1 1 dr) 2 \ l l 2 T{r) = th = const, n* Q = 1, respectively. Note that if n* Q and r were 
depending on two or three variables, a generalization of the routine GetRandomi) to 
more dimensions is possible. A less efficient possibility is to simulate r, n according 
to the element d 2 fdn and include the factor d 3 anQ/d 2 rdn in the residual weight in the 
step 6. Also note that the particle freeze-out coordinates calculated from relativistic 
hydrodynamics are distributed according to the element d 3 a tl u fJ- . 

4. The corresponding collective flow four- velocities u^(x) are calculated using, e.g., Eqs. 
(ED), dSHD or Eq. (J3HD- 

5. The particle three-momenta p*{sin 9 cos0, sin 9 sin 0, cos#} are then generated in the 
fluid element rest frames according to the probability f^ q (p°*;T, fii)p* 2 dp*dcos9*d(f)* 
by sampling uniformly distributed cos9* in the interval [—1,1], <fi* in the interval 
[0, 27r] and generating p* using a 100% efficient procedure similar to ROOT routine 
Get Random (). 

6. Next, the standard von Neumann rejection/acceptance procedure is used to account 
for the difference between the true probability W* i d 3 ad 3 p*/p°* (see Eqs. d20J), 
flUI) and the probability f- q (p°*;T,fXi)d 3 a^d 3 p* = f- q (p°*;T,fx i )n°*d 3 ad 3 f corre- 
sponding to the simulation steps 3-5. Thus the residual weight 

W*J 3 ad 3 p" , s 

W res = ^ (42) 

n°*p°*f^d 3 ad^ y } 

is calculated and the simulated particle four-coordinate and four-momentum are ac- 
cepted provided that this weight is larger than a test variable randomly simulated in 
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the interval [0, max(W/ es )]. Otherwise, the simulation returns to step 3. Note that for 
the freeze-out parametrizations considered in this paper, 



wr 



i 



n p 



(43) 



and the maximal weight max(W[ es ) can be calculated analytically. Particularly, in the 
Bjorken-like model and r) ma,x 3> 1, W? es is distributed in the interval [1 — tanh p™ ax , 1 + 
tanhp™ ax ]. The step 6 can be omitted for the Hubble-like model or for the Bjorken 
model without transverse flow (p u = 0) when W- es = 1. Generally, in the residual 
weight one should take into account the contribution of non-space-like sectors of the 
freeze-out hyper surf ace: 



1 - 



n*p* 

^ ip ! t" 



p n 



+ 9 



p n 



(44) 



7. Next, the hadron four- momentum p*^ is boosted to the fireball rest frame according 
to Eqs. (ESP- 

8. The two-body, three-body and many-body decays are simulated with the branching 
ratios calculated via ROOT utilities p|. A more correct kinetic evolution, taking into 
account not only resonance decays but also hadron elastic scattering, may be included 
with the help of the Boltzmann equation solver C++ code which was developed earlier 

It should be stressed that a high generation speed is achieved due to 100% generation 
efficiency of the freeze-out four-coordinates and four-momenta in steps 3-5 as well as due 
to a weak non- uniformity of the residual weight Wp s in the cases of practical interest. For 
example, in the Bjorken-like model, the increase of the maximal transverse flow rapidity from 



zero (Wl es = const) to p™ ax = 0.65 leads on 
speed. Compared, e.g., to THERMINATOR 
of magnitude faster. 



to a few percent decrease of the generation 
|, our generator appears more than one order 



VII. VALIDATION OF THE MC PROCEDURE 



In the Boltzmann approximation for the equilibrium distribution function ()14J) . i.e. re- 
taining only the first term in the expansion (fTT)|) . the transverse momentum (p t ) spectrum 
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in the Bjorken-like model takes the form 



dNi 1 /T f R fm t coshp u \ fp t sinhp u \ 

-g L t m t e Ml/ A77 / rdrKi I , (45) 



where ioC- 2 ) an d -^iX- 2 ) are the modified Bessel functions and m t = (rnj+p^) 1 ^ 2 is the particle 
transverse mass. 

To test our MC procedure, we compare in Fig. ^ the transverse momentum spectrum 
calculated according to Eq. (J45j) with the corresponding MC result for T = 0.165 GeV, 
R = 8 fm, mi = 0.14 GeV, Ar] = 10, \i = 0.0 GeV, r = 12 fm/c, p^ ax = 0.65 and 2.0. One 
may see that the analytical and the MC calculations practically coincide. 

To demonstrate the increasing influence of the residual weight with the increasing p™ ax , 
we also present in Fig. the MC results obtained without this weight. 



VIII. INPUT PARAMETERS AND EXAMPLE CALCULATIONS 

We present here the results of example MC calculations performed on the assumption of 
a common chemical and thermal freeze-out and compare them with the experimental data 
on central Au + Au collisions at RHIC. 



A. Model input parameters 

First, we summarize the input parameters which control the execution of our MC hadron 
generator in the case of Bjorken-like and Hubble-like parametrizations with a common ther- 
mal and chemical freeze-out: 

1. Number of events to generate. 

2. Thermodynamic parameters at chemical freeze-out: temperature (T) and chemical 
potentials per a unit charge (J1 b ,J1s,J1q). As an option, there is an additional pa- 
rameter 7 S < 1 taking into account the strangeness suppression according to partially 



equilibrated distribution 
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T, /*, 7 .) = — * (46) 

7 S exp dp* - Hi\/T) ± 1 
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TABLE I: The model parameters for central Au + Au collisions at \fs NN = 200 GeV. 



parameter 


Bj or ken-like 


Hubble-like 


T, GeV 


0.165 


0.165 


Jib, GeV 


0.028 


0.028 


Jls, GeV 


0.007 


0.007 


VQ, GeV 


-0.001 


-0.001 


7s 


1 (0.8) 


1 (0.8) 


r, fm/c 


6.1 


9.65 


R, fm 


10.0 


8.2 


'Jmax 


2 (3,5) 


2 (3,5) 


„max 
ru 


0.65 





where n| is the number of strange quarks and anti-quarks in a hadron i. Optionally, the 
parameter 7 S can be fixed usin g its phenomenological dependence on the temperature 
and baryon chemical potential 29]. 

3. Volume parameters: the freeze-out proper time (r) and firebal transverse radius (R). 



4. Maximal transverse flow rapidity (p™ ax ) for Bjorken-like parametrization 




5. Maximal space-time longitudinal rapidity (r/ max ) which determines the rapidity interval 
[ -I ?max) } ]mJ i n the collision center-of-mass system. To account for the violation of 
the boost invariance, we have included in the code an option corresponding to the 
substitution of the uniform distribution of the space-time longitudinal rapidity rj in 
the interval [— r] max , r] mauX ) by a Gaussian distribution exp(— rf /2/S.rf) with a width 
parameter At] (see, e.g., |3o|). 

The parameters used to model central Au+Au collisions at y/s NN = 200 GeV are given 
in Table |U 



B. Space-time distributions of the hadron emission points 

In figures 121 and EJ we show the distributions of the ir + emission transverse x-coordinate 
and time generated in the Bjorken-like and Hubble-like models with the parameters given 
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in Table HJ J)max = 2. Also shown are the contributions from the primary 7r + 's emitted 
directly from the freeze-out hypersurface and the contributions from 7r +, s from the decays 
of the most abundant resonances p, u, K*(892) and A. For primary pions, x < R and 
t < t < rcosh?7 max . The tails at \x\ > R and t > rcosh?7 max reflect the exponential law of 
the resonance decays. The longest tails in figures 121 and El are due to pions from u decays. 



C. Ratios of hadron abundances 



arge CTier^y 



It is well known that the particle abundances in heavy-ion collisions in a 
range can be reasonably well described within statistical models (see, e.g., 
based on the assumption that the produced hadronic matter reaches thermal and chemical 
equilibrium. This is demonstrated in tables |H] and 11111 for the particle number ratios near 
mid-rapidity in central Au +Au collisions at \/s ArAT = 130 and 200 GeV calculated in our 
MC model and the statistical model of Ref. [33j and compared with the RHIC experimental 



data. Being independent of volume and flow parameters, the particle number ratio allow one 
to fix the thermodynamic parameters. We have not tuned the latter here and simply used 
the same thermodynamic parameters as in Ref. despite there are noticeable differences 
in some particle number ratios calculated in the two models. These differences may be 
related to the different numbers of resonance states taken into account and uncertainties in 
the decay modes of high excited resonances. 



D. Pseudo-rapidity distributions 

In Fig. HJ we compare the PHOBOS data j^J on pseudo-rapidity spectrum of charged 
hadrons in central Au+Au collisions at y/s NN = 200 GeV with our MC results obtained 
within the Bjorken-like and Hubble-like models for different values of ?7 max . One may see 
that the data are compatible with the longitudinal boost invariance only in the mid-rapidity 
region in which the model is practically insensitive to ^ max - In the single freeze-out scenario, 
the data on particle numbers at mid-rapidity thus allows one to fix the effective volume 
V eS cx tR 2 . 
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TABLE II: Particle number ratios near mid-rapidity in central Au + Au collisions at y/s NN = 
130 GeV calculated with the thermodynamic parameters: T = 0.168 GeV, JIb = 0.041 GeV, 
Jl s = 0.010 GeV and JIq = -0.001 GeV. 



particle number ratios 


om MC 


statistical model [33] 


experiment 




0.98 


1.02 


1.00 ± 0.02 [34], 0.99 ± 0.02 [35J 


P/tt- 


0.06 


0.09 


0.08 ±0.01 [36] 


K-/K+ 


0.90 


0.92 


0.91 ± 0.09 [34], 0.93 ± 0.07 [37] 


K-/tt- 


0.22 


0.16 


0.15 ±0.02 [38] 


P/p 


0.61 


0.65 


0.60 ± 0.07 [34], 0.64 ± 0.08 [37] 


A/A 


0.69 


0.69 


0.71 ±0.04 [39] 




0.79 


0.77 


0.83 ± 0.06 [39] 


4>/K- 


0.17 


0.15 


0.13 ±0.03 [40J 


A/p 


0.48 


0.47 


0.49 ±0.03 [41], [42] 


ET fr- 


0.0086 


0.0072 


0.0088 ± 0.0020 [43J 



TABLE III: Particle number ratios near mid-rapidity in central Au + Au collisions at v^tvtv = 
200 GeV calculated with the thermodynamic parameters: T = 0.165 GeV, JIb = 0.028 GeV, 
/7s = 0.07 GeV, and JIq = -0.001 GeV. 



particle number ratios 



our MC 



experiment [45J 



7r /ir + 
K-/K J 

K-/TT- 
P/P 



0.98 
0.94 
0.21 
0.71 



0.984 ± 0.004 
0.933 ± 0.008 
0.162 ±0.001 
0.731 ±0.011 



E. Transverse momentum spectra 



In Fig. we compare the mid-rapidity PHENIX data 45] on n + , K + and proton p t 
spectra in Au±Au collisions at y/s NN = 200 GeV with our MC results obtained within the 
Bjorken-like and Hubble-like models. A good agreement between the models and the data 
may be seen for pions while for kaons and protons the models overestimate the spectra at 
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Pt < 1 GeV/c. For kaons, this discrepancy can be diminished with the help of the strangeness 
suppression parameter 7 S of 0.8 (see the right panel in Fig. |3J). The overestimated slope of 
the kaon and proton p t spectra can also be related with the oversimplified assumption of 
a common thermal and chemical freeze-out or insufficient number of the accounted heavy 
resonance states. 

The contribution of different resonances to the pion p t spectrum calculated in the Bjorken- 
like model is shown in Fig. |U1 

Note that in Hubble-like model, the transverse flow is determined by the volume param- 
eters R,t and so, at fixed thermodynamic parameters and the effective volume V e s oc tR 2 , 
the transverse momentum spectra allow one to fix both R and r. In the Bjorken-like model, 
there is more freedom since the transverse flow is also regulated by the parameter p™ ax . 
The choice of these parameters in Table |l] has been done to minimize the discrepancy of the 
simulated and measured correlation radii of identical pions (see below). 



F. Correlation functions 



It is well known that, due to the effects of quantum statistics (QS) and final state interac- 
tion (FSI), the momentum correlations of two or more particles at small relative momenta in 
their center-of-mass system are sensitive to the space-time characteristics of the production 
process on a level of fm = 10~ 15 m so serving as a correlation femtoscopy tool (see, for 
example, 5c|). 

The momentum correlations are usually studied with the help of correlation functions 
of two or more particles. Particularly, the two-particle correlation function CF(pi,p2) is 
defined as a ratio of the measured two-particle distribution to the reference one which is 
usually constructed by mixing the particles from different events of a given class, normalizing 
the correlation function to unity at sufficiently large relative momenta. 

Since our MC generator provides the information on particle four-momenta pi and four- 
coordinates Xi of the emission points, it can be used to calculate the correlation function 
with the help of the weight procedure, assigning a weight to a given particle combination 
accounting for the effects of QS and FSI. Here we will consider the correlation function of 
two identical pions neglecting their FSI, so the weight 

w = 1 + cos(g ■ Ax), (47) 
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where q = pi — P2 and Ax = x\ — x%. The CF is defined as a ratio of the weighted histogram 
of the pair kinematic variables to the unweighted one. 

Generally, the pair is characterized by six kinematic variables. In case of the azimuthal 
symmetry, there are five variables that are usually chosen as the three "out-side-long" com- 



ponents of the relative three-momentum vector [47|, |48J] q = (q out , q S ide, <?ion g ), half the pair 
transverse momentum k t and the pair rapidity or pseudo-rapidity. The out and side de- 
note the transverse, with respect to the reaction axis, components of the vector q; the out 
direction is parallel to the transverse component of the pair three-momentum. 

The corresponding correlation widths are usually parameterized in terms of the Gaussian 
correlation radii Ri, 

CF(p 1 ,p 2 ) = 1 + Aexp(-i?^ ut g^ ut - -R^de^Me - -Ri 2 ong <?i 2 ong - 2-R 2 ut;long g ut<?iong) (48) 

and their dependence on pair rapidity and transverse momentum is studied. The form of 
Eq. (plH|) assumes azimuthal symmetry of the production process j^. Generally, e.g., in 
case of the correlation analysis with respect to the reaction plane, all three cross terms q^qj 
contribute [3^. We choose as the reference frame the longitudinal co- moving system (LCMS) 
3|. In LCMS each pair is emitted transverse to the reaction axis so that the pair rapidity 
vanishes. The parameter A measures the correlation strength. For fully chaotic Gaussian 
source A = 1. Experimentally observed values of A < 1 are mainly due to contribution 
of very long-lived sources (rj, rj', K®, A, ...), the non-Gaussian shape of the correlation 
functions and particle misidentification. 

The correlation functions of two identical charged pions have been calculated within 
the Bjorken-like and Hubble-like models with the parameters given in Table HI ?7 max = 2, 
reasonably well describing single particle spectra in the mid-rapidity region. The three- 
dimensional correlation functions were fitted according to Eq. (pIKJ) in two kt intervals 0.1 < 
k t < 0.3 GeV/c and 0.3 < k t < 0.6 GeV/c. In Fig. [3 the fitted correlation radii and 
strength parameter are compared with those measured by STAR collaboration j^j]. One 
may see that the Bjorken-like model, adjusted to describe single particle spectra, describes 
also the decrease of the correlation radii with increasing k t but overestimates their values. 
The situation is even worth with the Hubble-like model which is more constraint than the 
Bjorken-like one and yields the longitudinal radius by a factor two larger. 

As for the overestimation of the correlation strength parameter A, it is likely related to 
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the neglected contribution of misidentified particles and pions from weak decays. Indeed, 
the new preliminary analysis of the STAR data with the improved particle identification [5l| 
yields the A parameter closer to the model results. 

We would like to emphasize that the high freeze-out temperature of 165 MeV and a fixed 
effective volume V e g oc tR 2 make it quite difficult to describe the correlation radii within 
the single freeze-out model. Thus a tuning of the longitudinal radius -Ri on g ~ TiT/mt) 1 ' 2 
requires a small proper time r, leading to too large values of R and i? S ide oc R. The concept 
of a later thermal freeze-out occurring at a smaller temperature T th < T ch and with no 
multiplicity constraint on the thermal effective volume (see section |H} can help to resolve 
this problem (see, e.g., 

To get a valuable information from the correlation data, one should consider more realistic 
models as compared with the simple Bjorken-like and Hubble-like ones (particularly, consider 
a more complex form of the freeze-out hypersurface taking into account particle emission 
from the surface of expanding system [20j] ) and study the problem of particle rescattering 
and resonance excitation after the chemical and/or thermal freeze-out ( only minor effect of 
elastic rescatterings on particle spectra and correlations is expected |25j). For the latter, 
our earlier developed C++ kinetic code 2^ can De coupled to the MC freeze-out generator. 



IX. CONCLUSIONS AND PERSPECTIVES 



We have developed a MC simulation procedure and the corresponding C++ code allow- 
ing for a fast but realistic description of multiple hadron production in central relativistic 
heavy ion collisions. A high generation speed and an easy control through input parame- 
ters make our MC generator code particularly useful for detector studies. As options, we 
have implemented two freeze-out scenarios with coinciding and with different chemical and 
thermal freeze-outs. Also implemented are various options of the freeze-out hypersurface 
parameterizations including those with non-space-like hypersurface sectors related to the 
emission from the surface of expanding system. The generator code is quite flexible and 
allows the user to add other scenarios and freeze-out surface parameterizations as well as 
additional hadron species in a simple manner. 

We have compared the RHIC experimental data with our MC generation results obtained 
within the single freeze-out scenario and Bjorken-like and Hubble-like freeze-out surface 
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parameterizations. While simplified, such a scenario nevertheless allows for a reasonable 
description of particle spectra. It however fails to describe the correlation functions of 
identical pions, overestimating the correlation radii. 

The RHIC data thus points to the need for a more complicated scenario likely including 
different chemical and thermal freeze-outs, a more complex form of the freeze-out hyper- 
surface (the use of numerical solution of the relativistic hydrodynamics is foreseen) and the 
account for kinetic evolution following the chemical and/or thermal freeze-out (for this, the 
MC generator can be coupled to our C++ kinetic code |25|). 

We plan to implement in the MC generator the impact parameter dependence of the 
freeze-out hypersurface and account for the anisotropic flow similar to Ref . 0, . In view of 
the importance of high-p t physics related to the partonic states created in ultra-relativistic 
heavy ion collisions, we also foresee the inclusion of mini-jet production 
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FIG. 1: The validation of the MC procedure for p^ = 0.65 (left panel) and 2.0 (right panel): the 
transverse momentum spectra (solid lines) calculated according to Eq. (|45[) and the corresponding 
MC results (black triangles). Also shown are the MC results obtained with a constant residual 
weight (black points). 




FIG. 2: The ir + emission transverse x-coordinate (left) and time (right) generated in the Bjorken- 
like model with the parameters given in Tabled imax = 2: all tt +, s (solid circles), direct 7r + 's (solid 
line), decay 7r + 's from p (squares), u> (open circles), K* (892) (up-triangles) and A (down-triangles). 
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FIG. 3: The same as in Fig. 21 for the Hubble-like parametrization. 




FIG. 4: The pseudo-rapidity (— lntan(0/2), 6 is the particle production angle) distributions of 
charged particles in central Au + Au collisions at y/s NN = 200 GeV from the PHOBOS experiment 
\a\ (solid circles) and the MC calculations within the Bj or ken-like (left panel) and Hubble-like 
(right panel) models. The model results corresponding to the space-time rapidity range parameter 
7?max = 5,3 and 2 are shown by solid, dashed and dotted lines respectively. 



25 




FIG. 5: The 7r + , K + and proton transverse momentum spectra at mid-rapidity y ~ in central 
Au + Au collisions at v^TViv = ^00 GeV from PHENIX experiment |45( (solid symbols) and the 
MC calculations within the Bjorken-like (dashed lines) and Hubble-like (solid lines) models. The 
right panel shows the model results obtained with the strangeness suppression parameter 7 S = 0.8. 
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p t (GeV/c) 

FIG. 6: The contributions to the tt + transverse momentum spectrum at mid-rapidity in central 
Au + Au collisions at y/s NN = 200 GeV calculated within the Bjorken-like model: all 7r + 's (solid 
circles), direct 7r + 's (stars), decay 7r + 's from p (squares), u> (open circles), K*(892) (up-triangles) 
and A (down-triangles). 
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FIG. 7: The ir^ir^ correlation radii and the suppression parameter A at mid-rapidity in central 
Au + Au collisions at y/s NN = 200 GeV from the STAR experiment £l| (open circles) and the MC 
calculations within the Bjorken-like model (up-triangles) in different intervals of the pair transverse 
momentum kt- 
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